#include <fstream>
#include <istream>
#include <sstream>
#include <vector>
#include <deque>
#include <string>
#include <iostream>



#define DEBUG_LEVEL 6
#include <debug_func.h>

#include <seq_pairs.h>
#include <blast_results.h>
#include <bcpp_read_seq.h>
#include <portable_timer.h>
#include <print_sequence.h>
#include <codeml.h>
#include <codeml_results.h>
#include <progress_indicator.h>
#define _MY_STL_MACROS
#include <stlhelper.h>
#include <ultoa.h>
#include <auto_restore_dir.h>
#include <print_error.h>

#include <boost/lexical_cast.hpp>
#include <boost/array.hpp>
#include <boost/foreach.hpp>
#include <boost/filesystem/path.hpp>
#include <boost/filesystem/operations.hpp>

#include "get_arg.h"


using std::vector;
using std::deque;
using std::ofstream;
using std::ifstream;
using std::istream;
using std::istringstream;
using std::string;
using std::binary_search;
using std::cerr;
using std::cin;
using std::cout;
using BioCPP::t_collect_sequences;
using BioCPP::t_collect_sequences_report_errors_to_strm;
using boost::filesystem::current_path;
using boost::filesystem::path;
using boost::lexical_cast;
using bcpp_align::t_align_results;

struct t_align_results_pair
{
	unsigned seq_id1;
	unsigned seq_id2;
	t_align_results blast_result;
	t_align_results_pair(){}
	t_align_results_pair(unsigned n1, unsigned n2, const t_align_results& blast_result_)
			:seq_id1(n1), seq_id2(n2), blast_result(blast_result_){}
};


//________________________________________________________________________________________

//		insert_gaps_into_sequence_pair

//________________________________________________________________________________________
void insert_gaps_into_sequence_pair(	const t_align_results&		pair_data,
										string&						sequence1,
										string&						sequence2,
										bool						cdna,
										const string&				pair_identifier)
{
	debug_func_cerr(5);
	unsigned codon = cdna ? 3 : 1;
	try
	{

		if (!pair_data.align_len)
			throw std::runtime_error("Zero length alignment!");


		if (!sequence1.length() || !sequence2.length())
			throw std::logic_error("Sequences missing.");

		if (cdna)
		{
			if (!expand_gaps_cdna(pair_data.gaps1, sequence1) )
				throw std::runtime_error("Alignment gaps do not match sequence1.");
			if (!expand_gaps_cdna(pair_data.gaps2, sequence2) )
				throw std::runtime_error("Alignment gaps do not match sequence2.");
		}
		else
		{
			if (!expand_gaps_prot(pair_data.gaps1, sequence1) )
				throw std::runtime_error("Alignment gaps do not match sequence1.");
			if (!expand_gaps_prot(pair_data.gaps2, sequence2) )
				throw std::runtime_error("Alignment gaps do not match sequence2.");

		}
		if (sequence1.length() != sequence2.length())
			throw std::runtime_error("Sequence lengths mismatch after inserting gaps.");

		if (pair_data.align_beg * codon > sequence1.length())
			throw std::runtime_error("Align start too large.");

		// remove non-aligned region
		sequence1.erase(0, pair_data.align_beg * codon);
		sequence2.erase(0, pair_data.align_beg * codon);

		if (pair_data.align_len * codon > sequence1.length())
			throw std::runtime_error("Align length too large.");

		sequence1.erase(pair_data.align_len * codon);
		sequence2.erase(pair_data.align_len * codon);
		return;

	}
	catch (std::exception& e)
	{
		throw std::runtime_error(pair_identifier + e.what());
	}
	catch (...)
	{
		throw std::runtime_error("Unknown gap insertion error");
	}

}




//________________________________________________________________________________________

//		do_parse_pairs

//________________________________________________________________________________________
struct do_parse_pairs
{
	const t_collect_sequences&		seq_db;
	vector<t_align_results_pair>& 		alignment_data;
	unsigned						name1_index;
	unsigned						name2_index;
	t_program_args&					args;
	std::string						comment_letters;
	bool							name_matched;
	do_parse_pairs(	const t_collect_sequences&			seq_db_,
					vector<t_align_results_pair>& 			alignment_data_,
					unsigned							name1_index_,
					unsigned							name2_index_,
					t_program_args&						args_,
					const std::string&					comment_letters_):
			seq_db(seq_db_), 
			alignment_data(alignment_data_),
			name1_index(name1_index_),
			name2_index(name2_index_),
			args(args_),
			comment_letters(comment_letters_),
			name_matched(name1_index == unsigned(-1) ? true : false){}
	void print_error(unsigned cnt_line, const string& prot_name) const
	{
		std_print_error(args.CERR(), "The sequence for \"" + prot_name +
								"\" on line #" + ultoa(cnt_line) +
								" could not be found");
	}
	void operator()(istream& is, unsigned cnt_line)
	{
		// get data for one protein pair
		string name1, name2;

		// comment if false
		string comment_line;
		if (!seq_pairs::parse_ids(is, name1, name2, comment_line, comment_letters,
								  seq_pairs::ePARSE_IDS_ONLY))
			return;


		
		if (args.kaks)
		{
			t_codeml_results ignore_kaksdata;
			is >> ignore_kaksdata;
		}
		t_align_results aligndata;
		is >> aligndata;

		if (!is)
			throw std::runtime_error("The alignment data is not in the correct format");
		
		// if names are already in sequence file, must be indexed
		unsigned index1, index2;
		if (!seq_db.get_index_of_name(name1,  index1))
		{
			print_error(cnt_line, name1);
			alignment_data.push_back(t_align_results_pair(unsigned(-1), unsigned(-1), aligndata));
			return;
		}
		if (!seq_db.get_index_of_name(name2,  index2))
		{
			print_error(cnt_line, name2);
			alignment_data.push_back(t_align_results_pair(unsigned(-1), unsigned(-1), aligndata));
			return;
		}

		alignment_data.push_back(t_align_results_pair(index1, index2, aligndata));

		// this  alignment matches names, put that on the to print list
		if ((index1 == name1_index && index2 == name2_index) ||
			(index2 == name1_index && index1 == name2_index))
		{
			args.pair_indices.push_back(alignment_data.size() - 1);
			name_matched = true;
			args.VERBOSE() << "\tAlignment on line #" << cnt_line << " matched\n";
		}
	}
};


		




void read_sequences(t_collect_sequences& seq_db, t_program_args& args)
{
	{
		ifstream f_cdna_data(args.file_seq.c_str());
		if (!f_cdna_data)
			throw std::runtime_error("Could not open " + args.file_seq +
										" to retrieve cDNA sequences");

		args.VERBOSE() << "\tReading ";
		if (args.cdna)
			args.VERBOSE() <<  "cDNA ";
		else
			args.VERBOSE() <<  "protein ";
		args.VERBOSE() <<  "sequences from [ " << args.file_seq << " ]...\n";
		t_progress_indicator dots(args.VERBOSE(), 1000);
		t_collect_sequences_report_errors_to_strm report_errors(args.CERR());
		seq_db.init_from_stream(f_cdna_data, args.regex, dots, report_errors);
	}
	args.VERBOSE() << "\tDone...\n";
	args.VERBOSE() << "\t" << seq_db.size() << " sequences\n";
}

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//		Main

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
int main(int argc, char* argv[])
{
	t_program_args args;

	try
	{

		debug_func_cerr(5);

		// don't use C streams
		std::ios::sync_with_stdio(false);


		// get command line arguments
		if (!process_command_line_options(argc, argv, args))
			return 1;


		//
		// read sequences
		//
		t_collect_sequences seq_db;
		read_sequences(seq_db, args);

		unsigned name1_index;
		unsigned name2_index;
		// do not search using the names if unspecified
		if (args.name1.length() == 0)
		{
			// will never match to -1
			name1_index =
			name2_index = unsigned(-1);
		}
		else
		{
			if (!seq_db.get_index_of_name(args.name1,  name1_index))
				throw std::runtime_error("The sequence for " + args.name1 + " could not be found.");
			if (!seq_db.get_index_of_name(args.name2,  name2_index))
				throw std::runtime_error("The sequence for " + args.name2 + " could not be found.");
		}

		//
		// read sequences alignment data
		// 		the indices of the specied names in args.name1 and name2
		// 		are added to pairs_indices as well
		//
		args.VERBOSE() << "\tFinding alignment data...\n";
		vector<t_align_results_pair> alignment_data;
		do_parse_pairs parse_pairs(seq_db,  alignment_data, name1_index, name2_index, args,
									args.comment_letters);
		seq_pairs::read_data(parse_pairs, "STDIN", args.VERBOSE());
		if (!parse_pairs.name_matched)
			throw std::runtime_error("No alignment data matching the identifiers \"" +
									args.name1 + "\" and \"" +
									args.name2 + "\" could be found.");


		//
		// align
		//
		if (args.pair_indices.size())
		{
			vector<unsigned>::const_iterator i = args.pair_indices.begin();
			vector<unsigned>::const_iterator end = args.pair_indices.end();
			for (;i != end; ++i)
			{
				// beyond range
				if (*i >= alignment_data.size())
					continue;

				// not in sequence file
				if (alignment_data[*i].seq_id1 == unsigned(-1) ||
					alignment_data[*i].seq_id2 == unsigned(-1))
					continue;

				string seq1 = seq_db.sequence(alignment_data[*i].seq_id1);
				string seq2 = seq_db.sequence(alignment_data[*i].seq_id2);
				string acc1 = seq_db.accession(alignment_data[*i].seq_id1);
				string acc2 = seq_db.accession(alignment_data[*i].seq_id2);
				insert_gaps_into_sequence_pair(alignment_data[*i].blast_result, seq1, seq2,
												args.cdna, "For " + acc1 + " vs " + acc2 + ".\n");


				if (args.wrap)
				{
					print_fasta(cout, acc1, seq1, 60);
					print_fasta(cout, acc2, seq2, 60);
				}
				else
					cout
						<< ">" << acc1 << "\n" << seq1 << "\n"
						<< ">" << acc2 << "\n" << seq2 << "\n";
			}
		}
		else
		{
			BOOST_FOREACH(t_align_results_pair& alignment, alignment_data) 
			{
				if (alignment.seq_id1 == unsigned(-1) ||
					alignment.seq_id2 == unsigned(-1))
					continue;
				string seq1 = seq_db.sequence( alignment.seq_id1);
				string seq2 = seq_db.sequence( alignment.seq_id2);
				string acc1 = seq_db.accession(alignment.seq_id1);
				string acc2 = seq_db.accession(alignment.seq_id2);
				insert_gaps_into_sequence_pair(alignment.blast_result, 
												seq1, seq2,
												args.cdna, 
											   "For " + acc1 + " vs " + acc2 + ".\n");

				if (args.wrap)
				{
					print_fasta(cout, acc1, seq1, 60);
					print_fasta(cout, acc2, seq2, 60);
				}
				else
					cout
						<< ">" << acc1 << "\n" << seq1 << "\n"
						<< ">" << acc2 << "\n" << seq2 << "\n";
				
			}
		}



	}
	catch (std::exception& e)
	{
		std_print_error(args.CERR(), e.what());
		return 1;
	}

	return 0;
}

